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We investigate the effect of hydrogen on the mobility of a screw dislocation in body-centered cubic (bee) iron 
using first-principles calculations, and show that the increase of screw dislocation velocity is expected in a limited 
temperature range. The interaction energy between a screw dislocation and hydrogen atoms is calculated for 
various hydrogen positions and dislocation configurations with careful estimations of the finite size effects, and 
the strongest binding energy of a hydrogen atom to the stable screw dislocation configuration is estimated as 
256 ± 32 nieV. These results are incorporated into a line tension model of a curved dislocation line to elucidate 
the effect of hydrogen on the dislocation migration process. Both the softening and hardening effect of hydrogen, 
caused by the reduction of kink nucleation enthalpy and kink trapping, respectively, are evaluated. A clear 
transition between softening and hardening behavior at the lower critical temperature is predicted, which is in 
qualitative agreement with the experimental observation. 

Keywords: First-principles calculation; Dislocations; Hydrogen Embrittlement; Hydrogen Enhanced Localized 
Plasticity 



1. Introduction 

Interaction between dislocations and solute 
atoms plays a key role in the solute hardening 
of metals, and in the case of body-centered cubic 
(bcc) metals sometimes solute softening occurs 
at low temperatures |l|2j . Solute atoms either 
pin dislocations and hinder their motion, or re- 
duce the Peierls barrier locally and promote kink 
nucleation and thus enhance dislocation motion. 
The precise effect of the solute atoms on plastic- 
ity depends on the core structure of the disloca- 
tion and binding energy landscape of the solute 
atom around the core, and a direct investigation 
of such atomistic scale properties necessarily re- 
quires first-principles calculations. In the present 
work we use the density functional theory (DFT) 
calculations to investigate the effect of hydrogen 
on the mobility of a dislocation in bcc Fe. 



Hydrogen is a special solute element because of 
its ubiquity in the environment and permeability 
into metals, resulting in a unique phenomenon 
known as hydrogen embrittlement (HE) in which 
fracture toughness of the material is reduced sig- 
nificantly when subjected to hydrogen-rich envi- 
ronment [5]. Various mechanisms of HE have 
been proposed, such as hydrogen enhanced de- 
cohesion [4 5 6J, suppression of dislocation emis- 
sion by hydrogen [7], and hydrogen enhanced lo- 
calized plasticity (HELP) [8]. Each mechanism 
qualitatively contributes to HE, and the degree 
of its contribution is expected to be sensitive to 
the strain rate, hydrogen density, and tempera- 
ture. Therefore, a quantitative estimate of the 
contribution of each mechanism is important to 
identify the dominant mechanism under realistic 
environmental and engineering conditions and to 
develop a material less susceptible to HE. 
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Softening of materials by hydrogen (H) solute 
atoms has been ascribed to a possible cause of HE 
in the HELP mechanism, based on the observa- 
tion of the reduced flow stress [^ and increased 
screw dislocation velocity [1^ when H atoms are 
introduced into bcc Fe, although hardening by H 
atoms also occurs depending on the subtle dif- 
ference in experimental settings [8]. In bcc Fe, 
H atoms concentrate on the stretched region un- 
der the tensile conditions [11I12| . If concentrated 
H atoms induce local slips, the local strain and 
dislocation density increase and further concen- 
trations of H atoms will occur, ultimately leading 
to the plastic instability and ductile fracture. 

Since the mobility of a dislocation in bcc met- 
als is mainly determined by that of the screw 
component owing to its large Peierls barrier and 
slow migration |13j . we investigate in the present 
work the interaction between a screw dislocation 
and an H atom. Although an edge dislocation 
strongly attracts H atoms by its long range hydro- 
static strain field and thus the effect of hydrogen 
on the mobility is expected to be much stronger 
compared to a screw dislocation, the mobility of 
an edge dislocation in bcc metals is very high and 
H atoms are likely to lower the mobility in most 
conditions [M]. On the other hand, the mobility 
of a screw dislocation in low temperatures is con- 
trolled by the rate of atomic scale kink nucleation, 
and a single H atom on a long screw dislocation 
segment can affect its nucleation rate. 

Fig. [T] shows core structures of a screw dis- 
location for several different core positions in 
bcc crystals, identified by the DFT calculations 
|15I16|17I18|I91 . The "easy core" configuration 
(EGG) shown in Fig. [T] (a) is the most stable 
configuration, while the "hard core" configura- 
tion (HGG) shown in Fig. [T] (b) is unstable or 
metastable and has a higher core energy owing to 
the large free volume inside the core. Movement 
of a screw dislocation in any direction requires the 
alternation of core structures between the EGG 
and HGG, resulting in the large Peierls barrier. 
Since H atoms are attracted to a free volume and 
lower the total energy, the trapping of an H atom 
to the HGG is expected to be stronger than that 
of the EGG, and the Peierls barrier is lowered by 
the H atom. 



In the case of a screw dislocation in bcc Fe, pre- 
vious DFT calculations have shown that the sad- 
dle point of a migration path between two adja- 
cent EGGs is close to, but different from the HGG, 
as shown in Fig. [T](c), and that the hard core and 
saddle point configuration (SPG) have nearly the 
same core energy 119!. Thus we assume that the 
saddle point moves toward the HGG when an H 
atom is in its core, and we calculate hydrogen 
binding energy for the EGG and HGG to investi- 
gate the effect of an H atom on the Peierls barrier. 

Although the effect of an H atom to lower the 
Peierls barrier seems obvious when an H atom 
is just ahead of the screw dislocation on the slip 
plane, a question arises whether or not the av- 
erage mobility is increased. Fig. [2] depicts the 
effects of an H atom on the migration process of 
a screw dislocation in bcc Fe. The migration is 
initiated by the thermal activation of a kink pair 
nucleus, followed by the movement of the kinks 
to the end of the straight screw dislocation seg- 
ment, or annihilation with other kinks. For the 
increased average mobility, at least one H atom 
must always be close to and ahead of a screw dis- 
location line to promote the kink pair nucleation, 
as shown in Fig. [5] (a). In addition, the H atom 
just behind the screw dislocation line can slow or 
stop the kink motion and decrease the dislocation 
mobility, as shown in Fig. [2] (b). In the present 
work, both of these competing effects are evalu- 
ated using a line tension model which is based 
on the DFT calculations, and conditions for the 
hydrogen softening are derived. The binding en- 
ergy between hydrogen and a screw dislocation 
have already been calculated by DFT in Ref . [20] , 
but the estimate in [30] does not include the zero 
point energy corrections and based on small num- 
ber of k-point samplings, which makes it difficult 
to compare the results with experiments. In the 
present work, we obtain more reliable estimates 
and compare the results with the experimental 
observations. 

The rest of this paper is organized as follows. 
We first describe the method used for the DFT 
calculations, including the details on the bound- 
ary conditions. In Section 3 the results of the 
DFT calculations are summarized. In Section 4 
the softening and hardening effect of hydrogen is 
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estimated by the line tension model, and the re- 
sults are compared with experiment. Concluding 
remarks are given in Section 6. 

2. Details of the DFT calculations 

The electronic structure calculations and the 
structure relaxations by force minimizations in 
the DFT steps are performed using the Vi- 
enna Ab-initio Simulation Package (VASP) pTl 
122] with the projector augmented wave method 
and ultrasoft pseudopotentials. The exchange 
correlation energy is calculated by the generalized 
gradient approximation with the Perdew-Burke- 
Ernzerhof function [53]. Spin-polarized calcula- 
tions are employed in all cases. The Methfessel- 
Paxton smearing method with 0.1-eV width is 
used. The cutoff energy for the plane-wave basis 
set is 350 eV, and the convergence of hydrogen so- 
lution energy with respect to the increasing cutoff 
is confirmed. Structural relaxation is terminated 
when the maximum force acting on the movable 
degrees of freedom becomes less than 10 meV/A. 

The reference configurations of a screw disloca- 
tion without an H atom is obtained using a modi- 
fied version of the fiexible boundary method [TS] . 
It consists of DFT relaxation of atoms in the core 
region and Green function relaxation of atoms 
outside the core region. In the DFT calculations, 
the system is divided into two concentric hexag- 
onal regions 1 and 2, and the atoms in region 1 
are relaxed while the atoms in region 2 are fixed. 
In the subsequent Green function relaxation, each 
atom's displacement relative to the linear elastic 
solution 24^ is calculated (denoted by Ui). Then 
a large number of atoms are added outside the 
region 2, whose positions are given by the linear 
elastic solution. Displacement field Ui of atoms in 
region 2 and outside are relaxed so that forces act- 
ing on these atoms calculated by the Hessian ma- 
trix of a perfect crystal, fi = J2iM -^ii('"j ~ ^0: 
becomes zero. The minimum image convention 
is applied for the difference of the displacement 
to account for the periodicity of the lattice. Dis- 
placements in the region 1, which is induced by 
the "core force" J25j, is fixed in this step. The 
DFT and Green function relaxation steps are re- 
peated until convergence. The number of atoms 



in region 1 and 2 are 48 and 99, respectively in 
the present work, as shown in Fig. [S] Through- 
out the present paper, Cartesian coordinates X, 
Y and Z which are parallel to (211), (Oil) and 
(111) , respectively (see Fig. [S]) are employed. 
The cell edge of the Z direction is equal to the 
Burgers vector, whose length is b. 

The two kinds of configuration ECC and 
HCC are used as reference configurations. For 
the hexagonal supercell, k-points are placed 
on a Gamma-centered mesh in the XY-plane 
to preserve the hexagonal symmetry and the 
Monkhorst-Pack k-point mesh is used in the Z 
direction. The numbers of k-points will be shown 
later for each calculation case. The lattice con- 
stant is set to qq = 2.833A, which is obtained by 
the DFT calculation of two Fe atoms with 20^ 
k-points. 

After the reference state is obtained, an H atom 
is placed at each tetrahedral site (t-site) near the 
core and the atomic configurations are relaxed for 
each case. The hydrogen solution energy is calcu- 
lated as Es ^ Ed+H -Ed- ^i/2/2, where Ed+n, 
Ed, and Eh2 denote energies of a dislocation with 
an H atom, a reference dislocation configuration 
and a hydrogen molecule, respectively. The bind- 
ing energy Ei, of a specific site is defined as a 
difference of Eg with respect to the bulk t-site so- 
lution energy, and defined as positive when Es is 
lower than that of the bulk t-site. 

The zero point energy (ZPE) correction to the 
solution energy is calculated from the Hessian ma- 
trix eigenvalues. We assume that Fe atoms are 
heavy enough compared to the H atom so that 
ZPE can be approximated by the motion of an 
H atom only. The corresponding Hessian matrix 
is calculated by displacing the H atom in each of 
the ±X, ±Y and ±Z directions by 0.015A and 
observing the force acting on the H atom. ZPE 
correction is calculated as 

^-^^E^Vfc^A^, (1) 

1=1 
where ki are the three eigenvalues of the 3x3 Hes- 
sian matrix, mn is the mass of an H atom and h 
is the Planck constant. We have confirmed that 
convergence of Ez with respect to the cut-off en- 
ergy, k-point mesh size and system size is fast and 
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errors are on the order of a few meV, which is neg- 
Ugible compared to other sources of errors. ZPE- 
corrected solution energy and binding energy are 
denoted by E^ and E^ , respectively. £'f is given 
by Sf = E, + E,~ E,H2/2, where E,h2 = 266 
meV is the ZPE of H2 molecule [26] . 

In the present work, the reference dislocation 
configuration is replicated and three layers are 
stacked in the Z direction to isolate the H atom 
from its mirror images. The number of atoms 
in this configuration is 341, which is beyond the 
limit of feasible computation in VASP. However, 
not all Fe atoms are necessary to calculate the 
hydrogen solution energy, since the displacement 
of the Fe atoms caused by the hydrogen solu- 
tion rapidly attenuates as the distance from the 
H atom increases. In the present work, Fe atoms 
in a hexagonal region that is centered at the H 
atom are clipped out from the reference configu- 
ration and used for the calculation of the binding 
energies, as shown in Fig. |4l Outer Fe atoms 
are fixed to the reference configuration, because 
these atoms are exerted by the strong artificial 
force caused either by a vacuum region or a do- 
main boundary [TS] . In the DFT calculations of a 
single dislocation, the domain boundary is gener- 
ally preferred over the vacuum boundary since the 
vacuum region induces a large amount of charge 
redistribution and a strong surface effect. How- 
ever, if the domain boundary condition is used, 
the k-point sampling in the XY-plane is required, 
while for the vacuum boundary only one k-point 
in the XY-plane is suffice since wave functions 
of any Bloch wave vector are allowed. In the 
present study, the domain boundary condition is 
employed and 3x3x8 k-points are used for the 
high-precision calculations, while for other cases 
1x1x8 k-points are used. The numerical errors 
caused by the k-point sampling for the 3x3x8 
and 1x1x8 cases are estimated as 5 meV and 
30 meV, respectively. 

Here, it is crucial to estimate the sufficient sys- 
tem size for the calculation of E^ with reasonably 
small finite-size errors. For this purpose, the fi- 
nite size effect on Eg is calculated using cubic bcc 
cell with 2L^ Fe atoms with L = 2, 3 and 4 with 
an H atom placed in the t-site. The Monkhorst- 
Pack k-point mesh of width 7r/12ao is used for 



each case. The solution energy without the re- 
laxation of Fe atoms is 441, 456 and 442 meV 
for L = 2, 3 and 4, respectively. As has been 
shown in [26] , Eg without relaxation quickly con- 
verges in the small size, so that the finite size ef- 
fect mainly comes from the relaxation. Fig. [5] (a) 
shows the energy change by the relaxation plot- 
ted against the inverse of the number of atoms 
under the constant-volume boundary conditions. 
Fig. [5] (a) also shows relaxation energy calcu- 
lated using the embedded atom method (EAM) 
potential which is fitted to various properties of 
hydrogen in the bcc Fe crystal calculated by DFT 
[27] . This EAM potential is fitted to energies 
without ZPE corrections and is suitable for the di- 
rect comparison with DFT calculations, although 
for the comparison with experiment, path integral 
molecular dynamics [27] or ZPE-corrected poten- 
tial '28' is required. One can see that, in both 
DFT and EAM cases, the finite-size effect of the 
relaxation energy is inversely proportional to the 
system size. 

Fig. [5] (b) shows the size dependence of the 
total pressure induced by the H atom solution, 
which is also inversely proportional to the sys- 
tem size (the residual pressure of the DFT case 
in the large volume limit is an artifact caused by a 
small error of the lattice constant of about 0.1%). 
This size dependence can be explained by assum- 
ing that the displacements of the Fe atoms caused 
by the H atom solution are inversely proportional 
to the square of the distance from the H atom. 
Es converges to the limit 188 meV as the system 
size increases while the cell volume is fixed, which 
is consistent with the DFT calculation with cell 
relaxation of Eg — 200 meV [3S]. From these 
results, a reasonable accuracy of 20 meV is ex- 
pected when the number of movable Fe atoms is 
about 30. In the present work, 144 Fe atoms in 
a three-layer hexagonal cell are used as a refer- 
ence configuration and only the 36 atoms in the 
inner hexagonal cell are relaxed, while the atoms 
in the outer buffer region are fixed. The width 
of the buffer region is large enough so that the 
maximum force induced by the domain boundary 
on the movable inner atoms is 0.03 eV/A for all 
cases. The finite size effect of this hexagonal cell 
is directly estimated from the perfect lattice t-site 
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solution energy using the same cell configuration 
shown in Fig. 21 The outer buffer atoms are fixed 
and 3x3x8 k-points are used, and the solution 
energy oi Eg = 215 meV is obtained. By com- 
paring it with the large volume limit 188 meV, 
the finite size error is estimated as 27 meV. To- 
gether with the error from the k-point sampling, 
the overall precision of Es is estimated as 32 meV 
and 57 meV for 3 x 3 x 8 and 1x1x8 k-points 
cases, respectively. The error for E^, is expected 
to be much smaller because of the cancellation of 
errors, but it is safe enough to employ the same 
error-bars for Ei,. 

As mentioned in Seel, the hydrogen binding 
energy of the HCC is expected to be larger than 
that of ECC. In DFT calculations, the hydrogen 
density inevitably becomes very high owing to the 
system size limitations and the ECC is expected 
to change its core structure to the HCC to gain 
more hydrogen binding energy at the expense of 
some core energy. However, to construct a model 
of the hydrogen-dislocation interaction, the esti- 
mation of the hydrogen binding energy in the di- 
lute hydrogen concentration limit (where the core 
structure is not affected by the hydrogen solu- 
tion) is required. To calculate the dilute hydro- 
gen concentration limit, the core position is fixed 
to the reference configuration at a position half 
the thickness of the system away in the Z direc- 
tion from the H atom. Since the core position 
is related to the Z displacements of three atoms 
surrounding the core [19], the core position at a 
certain layer can be fixed by prohibiting relax- 
ation in the Z direction of three atoms around 
the core. This way, the dislocation core remains 
in the reference configuration and the dilute con- 
centration limit of Es is obtained. 

3. Results 

Fig. HI shows the positions of binding sites 
obtained by the DFT calculations for ECC and 
HCC. The binding energies are summarized in 
Table [TJ Fig. |6] also shows minimum values of 
electron density along the screw dislocation line, 
which can be used to locate free volumes. As is 
expected, the strong binding sites are located in 
the regions of free volume indicated by the low 



electron density. The strongest binding site for 
ECC is E2, which is located around the triangle 
adjacent to the core, and the binding energy is 
£'f = 256 ± 32 meV. The binding site El has 
slightly smaller binding energy of 249 meV, and 
one can see that there are three broad potential 
basins that encompass the binding sites El and 
E2. The four binding sites (two El sites and two 
E2 sites) in a basin are close together and the 
energy barriers between them are expected to be 
very low, and it is more appropriate to regard 
this basin as a single binding site which we refer 
to as E1/E2 basin hereafter. The configuration of 
the Fe atoms around the E1/E2 basin is shown in 
Fig. [T] These Fe atoms form a slanted triangular 
prism that has a free volume much larger than 
that of a perfect crystal, and the binding sites El 
and E2 are on the midplane of the prism. 

It is noteworthy that this broad basin had been 
predicted by the EAM calculations [27 , and the 
barrier between the binding sites in the basin is 
on the order of 20 meV in the EAM calculations. 
The binding energy of E2 with ZPE correction 
had also been calculated by two recent EAM po- 
tentials as 290 meV for potential B in Ref. [35] 
and 290 meV at 300 K in Ref. [27], both in good 
agreement with the present result. 

The equilibrium hydrogen concentration Cb at 
a binding site with a binding energy E^ is calcu- 
lated from the McLean's equation as follows: 

CoCxp{E^/kBT)/3 



a = 



l + Coexp(^f/fcBr)/3' 



(2) 



where Co is the bulk hydrogen concentration, ks 
is the Boltzmann constant and T is the temper- 
ature. The factor 1/3 comes from the fact that 
there are three t-sites per Fe atom. Fig. [5| shows 
temperature dependence of Cb at the E1/E2 basin 
for the two cases Co — 0.1 atom ppm (typical 
value in industrial environment) and Co = 10 
atom ppm (typical value in charged samples in ex- 
periments). One can see that the binding energy 
is strong enough to concentrate H atoms to the 
screw dislocation at room temperature. When Cb 
is on the order of 0.5, the H-H repulsion must be 
taken into account and Eq. ([2]) is not valid, so 
the plot only shows region of Cb less than 0.4. 
The strongest binding site for HCC is HO which 
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is located at the center of the core, and its binding 
energy is 390 ±32 meV. This is 130 meV stronger 
than the ECC case, indicating that an H atom 
can lower the Peierls barrier about 130 meV. This 
is consistent with the DFT calculation using the 
nudged elastic band method [20] in which bind- 
ing energy difference between the ECC and SPC 
is reported to be about 100 meV. The estimate of 
the actual reduction will be shown in the next sec- 
tion. It is also implied that when sufficiently large 
numbers of H atoms are trapped at the E1/E2 
basin, the dislocation core structure may change 
to HCC to gain more hydrogen binding energy 
at the expense of the core energy. If the core 
structure changes to HCC when the three E1/E2 
basins are occupied by H atoms with the same 
concentration Cb, H atoms at the three basins 
end up in the HO site and two H2 sites. The core 
energy increase by 40 meV per b [19] , and the to- 
tal hydrogen trap energy changes by an amount 
of C6(390 - 256) + 26*6(189 - 256) = meV per 
b on the average. Thus the ECC will remain un- 
changed regardless of the hydrogen concentration 
at the binding sites. 

It is clear that the effect of hydrogen on the 
migration of screw dislocation is strongest when 
one H atom is on the slip plane of the dislocation, 
as shown in Fig. [9l The relative position of the 
H atom with respect to the core changes as E8, 
H2, E2, HO, E2, H2 and E8 when the core moves 
from left to right, that is, the core approaches and 
leaves the H atom. Fig. [10] shows the hydrogen 
binding energies of these sites plotted against the 
distance r between the H atom and dislocation 
core. The DFT data is well fitted by a Lorentzian 
function 



390 

"""^^^ = l + 2(r/ro)^ 



meV, 



(3) 



where tq — \/6ao/3 is the distance between two 
neighboring ECC positions. Note that this func- 
tion is purely empirical one. Fig. [10] also shows 
prediction from the linear elasticity theory [TT] . 
As is expected, the binding energy at the core re- 
gion is much larger than that of the linear elastic 
interaction. 



4. Line tension model 

For the purpose of estimating the effect of H 
atoms on the kink pair nucleation enthalpy and 
kink nucleation rate, the interaction energy be- 
tween the H atom and the screw dislocation ob- 
tained in the previous section is incorporated into 
the line tension model of a dislocation line |29|30j , 
which is expressed as an enthalpy of a curved 
screw dislocation configuration specified by the 
core positions Pj at each atomic layer j of thick- 
ness b as follows: 



j 

+j2Ep{Pj) + {i<T*b)xr}-Pj 

~Y.Eh{\P,-PI^\), (4) 



where K — 0.866 eV/A^ is a constant derived 
from the Hessian matrix of ECC calculated by 
DFT [IH], Pj is a two-dimensional vector whose 
components are the X and Y coordinate of the 
dislocation core position, Ep is the Peierls bar- 
rier per Burgers vector 6, Eh is the interaction 
energy between the dislocation line and the H 
atoms and P^ is the position of the fcth H atom 
in the XY-plane. The third term is the contri- 
bution from the external stress, where cr * 6 is 
a tensor- vector product of the stress and Burg- 
ers vector and I = Pj — Pj-i. The inclination 
of the dislocation line given by \Pj — Pj+i\/b is 
at most 1/30 [Slj, and we assume that the inter- 
action between the dislocation and the H atoms 
is described well by the hydrogen binding energy 
of a straight screw dislocation. Thus the inter- 
action term Eh in Eq. ([3]) is calculated only 
between an H atom and a representative dislo- 
cation segment which is closest to the H atom. 
In addition, we assume that the position of the 
H atom, which is initially placed at the E1/E2 
basin, remains virtually unchanged in the kink 
nucleation/migration process, since the binding 
in this site is strong throughout the process. The 
two dimensional Peierls energy Ep has been es- 
timated by the DFT calculations and fitted to a 
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plane- wave expansion as follows |19) : 



Ep{P^,Py) = ^{Ci/e(x,) + C2/o(x„) 



Q = l 

+ C5 [fe{xi -X2) + fe{x2 - X3) + fe{x3 - Xi)] , (5) 

where 

fe{x) = -(1 -cos27rx), 
foix) = -sin27rx, 

xi = APy, 

X2 = A{-Py + V3P^)/2, 
X3 = A{-Py-^P^)/2, 
A = V2/ao. 

The coefficients are Ci = 21.82, C2 = -14.51, 
C3 = 2.59, C4 = -2.72, and C5 = -2.89 meV. 

The kink nucleation enthalpy with H {Ekh) 
and without H {Ek) for a given shear stress ap- 
plied in the [111] (110) direction is calculated us- 
ing the string method fS? , for the dislocation mi- 
gration path E2-H0-E2. Fig. [TT] shows the stress 
dependence of Ekh and Ek, together with the 
experimental data of Ek from Ref. [33] . The re- 
duction of the enthalpy, AEk — Ek — Ekh, is 
about 110 meV for all the applied stress cases. 

The discrepancy in the stress dependence of 
Ek between the atomistic models and the ex- 
periments is a common problem for bcc metals 
[34|35|36j . and possible sources of this discrep- 
ancy are ascribed to several mechanisms, such as 
the pile-up effects [37] and the ZPE correction 
[55] . Nevertheless, since AEk only weakly de- 
pends on the stress and is mainly determined by 
the hydrogen binding energy, we assume that the 
estimate of AEk is valid for all stress regions. 

Fig. [12] (a) shows the kink nucle- 

ation/migration process in the presence of an 
H atom for the 400 MPa case. The dislocation 
line intersects with the H atom at the highest en- 
thalpy configuration shown by the bold line and 
lowers the enthalpy barrier. Fig. [12] (b) shows the 
highest enthalpy configuration seen from the Z 
direction. Contrary to our early expectation, the 



migration path of the dislocation is only slightly 
attracted to the HCC position from the SPC. 

At low temperature regime, the screw disloca- 
tion velocity is determined by the kink nucleation 
rate given by the following Arrhenius law: 



DdNdexpi-EK/keT), 



(6) 



where Dd is a prefactor and Nd is the length of 
the dislocation in unit of b. From the experimen- 
tal observation in Ref. [39], the nucleation rate 
of a dislocation of length Ndb = 2/im at 300 K 
and under the applied shear stress of about 33 
MPa is 81 s~^. The experimental value of Ek for 
the shear stress of 33 MPa is estimated from Fig. 
[TT] as 595 meV, and the prefactor in Eq. ([6]) is 
estimated as Dd = 0.99 x 10^ s"^ 

When the E1/E2 basins ahead of the disloca- 
tion in the slip direction are occupied by H atoms 
with the concentration Cf,, the nucleation rate is 
enhanced by a factor l+CbWk{exp{AEK /kBT) — 
1}, where Wk ^ 10 is the width of the kink nu- 
cleus in unit of b. Thus, for significant enhance- 
ment Cb must be on the order of 

r^ = I (7\ 

^' Wk{eM^EKlkBT)-l} ^'' 

or greater than this value. As the temperature 
is raised, C^ increases and Cb decreases. There 
is an upper critical temperature Tjj above which 
C^ > Cb is not satisfied. For the bulk hydrogen 
concentration of 0.1 and 10 appm, Tu is 280 and 
400 K, respectively. 

For the steady enhancement of the kink nu- 
cleation, the timescale of the hydrogen diffusion 
should be much shorter than that of the dislo- 
cation migration. After a dislocation migration 
event, the hydrogen concentration at the E1/E2 
basin ahead of the dislocation (referred to as 
"promotion site" hereafter) is much smaller than 
C^. The enhanced nucleation does not occur 
until the H atoms redistribute and the hydro- 
gen concentration at the promotion site becomes 
comparable with C^ . The jump rate of the H 
atom at the bulk has been calculated by DFT as 
Dnexpi-E'g/kBT), with Dh = 5.1 x lO^^ g-^ 
and E^ = 88 meV [T^, in reasonable agreement 
with the experiments [26| . Since the prefactor 
Dh is orders of magnitude greater than Dd, the 



The effect of Hydrogen atom on the Screw Dislocation Mobility in BCC Iron: A First- Principles StudyS 



timescale of the hydrogen redistribution is negh- 
gible if ExH is greater than E'^. From the exper- 
imental values of E^ in Fig. [TT] and the value of 
AEk = 110 meV obtained in the present work, 
the condition Ekh > Eh sets an upper critical 
shear stress aui = 190 MPa for the enhanced 
screw dislocation mobility, above which the hy- 
drogen diffusion cannot catch up with the dislo- 
cation motion [5]. We suppose that, below aui, 
the hydrogen diffusion is fast enough compared 
to screw dislocation motion so that the hydrogen 
density at the E1/E2 basin is approximated well 
by Ch in Eq. ^. For more precise evaluation, a 
kinetic Monte Carlo simulation of combined dislo- 
cation motion and hydrogen diffusion is expected. 

Next, we investigate the kink trapping effect of 
the H atom depicted in Fig. [131 Suppose that 
an H atom is at an E1/E2 basin behind the dis- 
location line, and a kink is moving from left to 
right. After the kink sweeps past the H atom, 
the relative position of H atom to the dislocation 
core changes from E1/E2 basin to E8, and its so- 
lution energy increases by 179 ± 57 meV. This 
kink trapping barrier, denoted by Et, is reduced 
as the applied shear stress a increases because 
the enthalpy as a function of the kink position z 
has a term —zabh, where h = v6ao/3 is the kink 
height, and there is a critical stress Cc at which 
the barrier vanishes, as shown in Fig. [T3](b). The 
values of Et for various applied shear stress are 
calculated using the string method, and are plot- 
ted against the shear stress in Fig. [141 together 
with the plots of Ek obtained by the DFT calcu- 
lations and experiments. The estimated value of 
ac is 265 MPa. 

Below CTc, the total time for a screw dislo- 
cation line to move to the adjacent Peierls en- 
ergy dip increases by the time required to escape 
from the traps of each H atom behind the dis- 
location. We assume that the average time of 
a single escape event follows the Arrhenius law 
with the same prefactor Dd of the kink nucle- 
ation as exp{Et/ ksT) / DdWk ■ The total escape 
time is then given by CbNdexp{Et/kBT)/DdWk- 
When this time is comparable to the aver- 
age nucleation time without hydrogen given by 
exp{EK /kBT)/NdDd, the effect of the enhanced 
nucleation rate is completely negated. Thus a 



condition 

CbNdexp{Et/kBT)_ exp{EK/kBT) 



DdWk 



NdDd 



(8) 



must be satisfied, which is simplified to Ex—Et < 
ksT \og{CbNj /Wk) ■ If we use the experimental 
values of Ek , the difference Ek — Et decreases as 
the shear stress is increased, and there is a second 
upper critical stress au2 above which Eq. ([5]) is 
not satisfied. 

When the hydrogen concentration at the 
E1/E2 basin is Cb, the enthalpy increase by the 
hydrogen de-trapping per unit length of kink mo- 
tion is CbE°/b, where E° = 179 meV. If this is 
larger than the enthalpy gain —abh, the kink mo- 
tion becomes impossible. Thus there is a lower 
critical stress ctl — CbEt/hb^ below which the 
screw dislocation cannot move. 

In total, there are four conditions which must 
be satisfied for the enhanced screw dislocation 
mobility of a screw dislocation, summarized as 
follows: 

1. The temperature must be below Tjj, above 
which hydrogen concentration at the pro- 
motion site is insuflficient for the en- 
hanced kink nucleation. Tu increases as 
Co increases. The boundary is given by 
C^{Tu) = Cb{Co,Tu). 

2. The applied shear stress must be greater 
than (Ti, below which a screw disloca- 
tion cannot move owing to the kink trap- 
ping by hydrogen, a^ increases as Cq in- 
creases. The boundary is given by ctl = 
CbiCQ,T)E?/hb\ 

3. The applied shear stress must be lower 
than (7(71 , above which dislocation motion is 
too fast and hydrogen concentration cannot 
catch up with it. am is independent of Cq. 
The boundary is given by Ekh{o'ui ) = E^^. 

4. The applied shear stress must be lower 
than (7(72, above which the time required 
for a kink to escape from H atoms behind 
the dislocation exceeds the average disloca- 
tion migration time without the H atom, 
cr [72 decreases as Cq increases or the length 
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of the screw dislocation increases. The 
boundary is given by £^(0112) — Et{cru2) = 
kBT\og{Ck{Co,T)NJlWk). 



Fig. [15] shows these conditions in the stress- 
temperature diagram for the two cases (a) Cq — 
0.1 appm and (b) Co = 10 appm, together with 
the experimentally observed yield stress compiled 
in Ref. [39]. In each figure, plots of au2 for 
two typical dislocation length, 2.0 and 0.2/im are 
shown. For the sake of comparison with the ex- 
periments, values of Ek are taken from the ex- 
periment to evaluate the four conditions. If we 
use the DFT values instead, the boundary of aui 
and au2 shifts to the right side. As can be seen 
in Fig. 1151 (t^ abruptly increases at tempera- 
ture where Cf, becomes on the order of 0.1, and 
exceeds the yield stress. Thus it is more appro- 
priate to regard the boundary marked by ul as a 
lower critical temperature. 

This result qualitatively agrees with the exper- 
imental observation of H-induced flow stress re- 
duction in the macroscopic samples by Matsui 
and Kitamura [9 . In their experiments, the bulk 
hydrogen concentration was estimated to be (at 
least) 15 — 30 appm. Reduction of the flow stress 
by the hydrogen charging was observed, and the 
reduction was more significant in lower tempera- 
tures (as is the general case of the solution soft- 
ening), whereas below 190 K, hardening by the 
hydrogen charging was observed. The maximum 
flow stress in their experiment was 200 MPa, thus 
the shear stress in any slip plane does not ex- 
ceed 100 MPa because the Schmid factor is al- 
ways less than 1/2. The transition at 190 K can 
be attributed to the crossing of 0-^, although for 
Co ~ 30 appm the crossing occurs at the much 
higher temperature in our model. The quantita- 
tive disagreement of the transition temperature 
can be attributed to either the numerical uncer- 
tainty of the binding energy, or to our approxima- 
tion of the hydrogen concentration at the binding 
site to the equilibrium concentration. The actual 
hydrogen concentration behind the moving dislo- 
cation must be much smaller, and in that case the 
transition temperature decreases. 



5. Summary and Conclusion 

Using DFT calculations with careful estima- 
tions of the finite size effects, we have calculated 
the binding energies of an H atom at various bind- 
ing sites around a screw dislocation of bcc Fe for 
two kinds of core configurations. The strongest 
binding energy to the stable and unstable core 
configurations were estimated as 256 ± 32 and 
390 ± 32 meV, respectively. Experimentally ob- 
served hydrogen binding energy to the screw dis- 
locations is currently not available, because in 
the thermal desorption spectroscopy, the desorp- 
tion peak for screw dislocations is insignificant 
compared to that of edge dislocations and grain 
boundaries. 

The interaction between an H atom and a screw 
dislocation was incorporated into a line tension 
model of a curved dislocation line. Using this 
model, the reduction of kink nucleation enthalpy 
by hydrogen was estimated to be 110 meV. The 
softening effect of H atoms by promoting the kink 
nucleation and the hardening effect by trapping 
the kink movement were both evaluated, and four 
conditions for the overall softening were derived. 
These conditions consist of two upper critical 
stresses and upper/lower critical temperatures. 
The upper critical stresses were found to be far 
above the yield stress, and only relevant for very 
high strain rate cases. The temperature range for 
the softening bounded by the upper/lower crit- 
ical temperatures was roughly estimated to be 
200 — 300 K for the bulk hydrogen concentration 
of 0.1 appm case and 300 — 400 K for 10 appm 
case. For the precise estimate of the temperature 
range, a kinetic Monte Carlo simulation which 
incorporates both the dislocation motion and the 
hydrogen diffusion is expected. A clear transition 
between softening and hardening behavior at the 
lower critical temperature is predicted, which is 
in qualitative agreement with the experimental 
observation f5^. 

While the dislocation migration behavior in fer- 
ritic steels is much different compared to the pure 
iron case owing to the presence of carbon and 
other solute atoms, the low temperature harden- 
ing caused by the dense hydrogen segregation to 
the screw dislocation is expected to be indepen- 
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dent of other solute atoms. It is expected that 
properties concerning the softening to hardening 
transition at low temperatures are observed for 
various hydrogen concentrations, strain rates and 
material purity. 

REFERENCES 

1. E. Pink and R. J. Arsenault, Prog. Mater. Sci. 
24, 1 (1979). 

2. D. R. Trinkle and C.Woodward, Science 9, 
1665 (2005). 

3. R. P. Gangloff, Critical Issues in Hydrogen 
Assisted Cracking of Structural Alloys, Else- 
vier Science, Oxford (2005). 

4. J. R. Rice and J. S. Wang, Mater. Sci. Eng. A 
107, 23 (1989). 

5. M. Yamaguchi, J. Kameda, K. Ebihara, 
M. Itakura and H. Kaburaki, Philos. Mag. 92, 
1349 (2012). 

6. J. J. Rimoh and M. Ortiz, Philos. Mag. 90, 
2939 (2010). 

7. J. Song and W. A. Curtin, Nat. Mater., 
doi:10.1038/nmat3479 (2012). 

8. H. K. Birnbaum and P. Sofronis, 
Mater.Sci.Eng. A 176, 191 (1994). 

9. H. Matsui and H. Kitamura, Mater.Sci.Eng. 
A 40, 207 (1979). 

10. T. Tabata and H. K. Birnbaum, Scripta 
Metal. 17, 947 (1983). 

11. S. Wang, K. Takahashi, N. Hashimoto, 
S. Isobe and S. Ohnuki, Scripta Mater. 68, 
249 (2013). 

12. A. Ramasubramaniam, M. Itakura, M. Or- 
tiz and E. A. Carter, J. Mater. Res. 23, 2757 
(2008). 

13. F. V. Vitek, Prog.Mater.Sci. 56, 557 (2011). 

14. S. Taketomi, R. Matsumoto and N. Miyazaki, 
J.Mater.Res. 26, 1269 (2011). 

15. C. Woodward and S. I. Rao, Phys. Rev. Lett. 
88, 216402 (2002). 

16. S. Frederiksen and K. Jacobsen, Philos. Mag. 
83, 365 (2003). 

17. E. Clouet, L. Ventelon and F. Willaime, Phys. 
Rev. Lett. 102, 055502 (2009). 

18. L. Ventelon and F. WiUaime, Philos. Mag. 90, 
1063 (2010). 

19. M. Itakura, H. Kaburaki and M. Yamaguchi, 



Acta Mater. 60, 3698 (2012). 

20. Y. Zhao and G.Lu, Model. Simul. Mater. Sci. 
Eng. 19, 065004 (2011). 

21. G. Kresse and J. Hafner, Phys. Rev. B 47, 
558 (1993). 

22. G. Kresse and J. Furthmiiller, Phys. Rev. B 
54, 11169 (1996). 

23. J. P. Perdew, K. Burke and M. Ernzerhof, 
Phys. Rev. Lett. 77, 3865 (1996). 

24. The anisotropic and isotropic linear elastic so- 
lutions of the displacement field of the screw 
dislocation are identical. See Ref. [19j . 

25. E. Clouet, Phys. Rev. B 84, 224111 (2011) 
E. Clouet, L. Ventelon and F. Willaime, 
Phys. Rev. B 84, 224107 (2011). 

26. D. E. Jiang and E. A. Carter, Phys. Rev. B 
70, 064102 (2004). 

27. H. Kimizuka and S. Ogata, Phys. Rev. B 84, 
024116 (2011). 

28. A. Ramasubramaniam, M. Itakura and 
E. A. Carter, Phys. Rev. B 79, 174101 (2009). 

29. D. Rodney and L. Proville, Phys. Rev. B 79, 
094108 (2009). 

30. K. Edagawa, T. Suzuki and S. Takeuchi, 
Phys. Rev. B 55, 6180 (1997). 

31. L. Ventelon, F. Willaime and P. Leyronnas, 
J. Nucl. Mater. 286-388, 26 (2009). 

32. W. E, W. Ren, and E. Vanden-Eijnden, Phys. 
Rev. B 66, 052301 (2002). 

33. W. A. Spitzig, Acta Metal. 18, 1275 (1970). 

34. C. Domain and G. Monnet, Phys. Rev. Lett. 
95, 215506 (2005). 

35. D. E. SegaU, A. Strachan, W. A. Goddard 
HI, S. Ismail-Beigi, and T. Arias, Phys. Rev. 
B 68, 014104 (2003). 

36. R. Groger, A. G. Bailey and V. Vitek, Acta 
Mater. 56, 5401 (2008). 

37. R. Groger and V. Vitek, Philos. Mag. Lett. 
87, 113 (2007). 

38. L. Proville, D. Rodney and M. C. Marinica, 
Nat. Mater. 11, 845 (2012). 

39. D. Caillard, Acta Mater. 58, 3493 (2010); 
ibid., 3504 (2010). 



The effect of Hydrogen atom on tlie Screw Dislocation Mobility in BCC Iron: A First-Principles Studyll 



Table 1 

Number of k-points Nk, solution energy E, 
to the binding energy AE^, total binding ener, 
binding site shown in Fig. [51 The "00" case is 
are in meV. 

Site Nk Es Eb Ez AEz 



binding energy Ef,^ ZPE correction E^, ZPE correction 

■gy, and estimated numerical and finite-size error for each 

the reference values of perfect crystal t-site. All energies 



E^ 



E,' 



Err 



00 


1 


X 


1 


X 


8 


241 





238 





346 





±57 




3 


X 


3 


X 


8 


215 





238 





320 





±32 


EO 


1 


X 


1 


X 


8 


156 


85 


224 


14 


247 


99 


±57 


El 


3 


X 


3 


X 


8 


22 


193 


182 


56 


71 


249 


±32 


E2 


3 


X 


3 


X 


8 


30 


185 


167 


71 


64 


256 


±32 


E3 


1 


X 


1 


X 


8 


103 


138 


175 


63 


145 


201 


±57 


E4 


1 


X 


1 


X 


8 


119 


122 


176 


62 


162 


184 


±57 


E5 


1 


X 


1 


X 


8 


211 


30 


245 


-7 


323 


23 


±57 


E6 


1 


X 


1 


X 


8 


231 


10 


248 


-10 


346 





±57 


E7 


1 


X 


1 


X 


8 


247 


-6 


234 


4 


348 


-2 


±57 


E8 


1 


X 


1 


X 


8 


195 


46 


207 


31 


269 


77 


±57 


E9 


1 


X 


1 


X 


8 


236 


5 


225 


13 


328 


18 


±57 


HO 


3 


X 


3 


X 


8 


-71 


286 


134 


104 


-70 


390 


±32 


HI 


3 


X 


3 


X 


8 


-56 


271 


187 


51 


-2 


322 


±32 


H2 


1 


X 


1 


X 


8 


116 


125 


174 


64 


157 


189 


±57 


H3 


1 


X 


1 


X 


8 


183 


58 


249 


-11 


299 


47 


±57 


H4 


1 


X 


1 


X 


8 


218 


23 


229 


9 


314 


32 


±57 


H5 


1 


X 


1 


X 


8 


236 


5 


222 


16 


325 


21 


±57 


H6 


1 


X 


1 


X 


8 


247 


-6 


227 


11 


341 


5 


±57 





Slip plane 



Slip plane 



Figure 1. Atomic structures of a screw dislocation 
in the (a) easy core, (b) hard core, and (c) migra- 
tion saddle point configurations, shown by gray 
spheres. White spheres are the atom positions in 
the perfect bcc crystal. 



Figure 2. Schematic diagram of the effects of 
an H atom on the migration process of a screw 
dislocation by (a) lowering the Peierls barrier and 
(b) decreasing the kink velocity. See the main 
text for details. 
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Figure 3. Two regions 1 and 2 used for the flexible 
boundary method in DFT, shown by white and 
black circles, respectively. The arrows show the 
difl'erentiated displacement of the easy core con- 
figuration. A Cartesian coordinate used through- 
out the present paper is also shown. 
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Figure 4. Calculation cell of the reference dis- 
location configuration (left) and dislocation+H 
configuration (right). Fe atoms close to the hy- 
drogen trap site is clipped out from the reference 
configuration and replicated in the Z direction 
three times, and H atom is placed at the trap site 
and structural relaxation is carried out. The cross 
symbol, white and black circles denote core posi- 
tion, movable and fixed Fe atoms, respectively. 
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Figure 5. Size dependence of (a) relaxation en- 
ergy and (b) pressure in the system of a bcc Fe 
perfect crystal with an H atom at the t-site under 
the constant-volume boundary conditions. The 
result of EAM potential [37] calculations is also 
shown. 
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Figure 9. A process of screw dislocation inter- 
action with H atom when the dislocation moves 
from left to right and the H atom is on the same 
slip plane. The labels denote the relative position 
of immobile H atom with respect to each disloca- 
tion position. 



Figure 7. Configuration of Fe atoms adjacent to 
the easy core. The cross marks the core position. 
The binding site El and E2 are on the midplane 
of the slanted triangular prism shown by the dark 
shaded triangle. 
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Figure 8. H concentration at the binding site 
close to the easy core screw dislocation, calcu- 
lated by McLean's equation with binding energy 
256 ± 32 meV for two cases of bulk hydrogen con- 
centrations. The dashed and dotted lines show 
the range of numerical uncertainty. 
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Figure 10. Hydrogen binding energy plotted 
against the distance between an H atom and dis- 
location core. A fitting curve to the DFT data is 
also shown by the solid line. The dashed line is 
the binding energy derived by the linear elasticity 
theory [llj . 



The effect of Hydrogen atom on tlie Screw Dislocation Mobility in BCC Iron: A First-Principles StudylA 



800 
700 
600 
500 
400 
300 
200 
100 







Ek-EkH 

Kperiment □ 




r-'Bi.E 


Y'v 


^\^ 






v/Sv,^ 






^^^i^ ^^\ 






V r "'---: 


^^^3^^ 







200 400 600 800 

Shear Stress (MPa) 



1000 



Figure 11. Kink nucleation enthalpy with H 
{Ekh) and without H {Ek) for the shear stress 
appUed in the [111] (110) direction. Estimated 
values of Ek from the experiment (Ref. [33]) are 
also shown with a linear fitting line. 




Figure 12. (a) Kink nucleation and migration 
process for 400 MPa case in the presence of an 
H atom shown by the white circle. The bold line 
represents the highest enthalpy configuration, (b) 
The highest enthalpy configuration seen from the 
Z direction. The positions of easy core, hard core 
and the saddle point of the Peierls barrier are also 
shown. 
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Figure 13. Schematic of kink trapping effect by 
the H atom, (a) As a kink moves from left to right 
and passes the H atom behind the dislocation line, 
the binding between the dislocation and the H 
atom is weakened, (b) The enthalpy increases by 
the amount of weakened binding energy. As the 
shear stress is applied, a slope proportional to the 
stress is added to the enthalpy curve and at some 
critical stress CTc the barrier vanishes. 
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Figure 14. Kink trap enthalpy Et plotted against 
the shear stress. Plots of the kink nucleation en- 
thalpy Ek is also shown. 
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Figure 15. A region in the temperature-stress 
diagram where the increase of screw dislocation 
velocity by hydrogen solution is possible for the 
bulk hydrogen concentration of (a) 0.1 appm and 
(b) 10 appm cases. The region is bounded by 
the upper critical temperature Tu, lower critical 
stress (Tl, the yield stress, and two kinds of upper 
critical stresses ajji and au2- <7U2 depends on 
the dislocation length, and two cases of typical 
dislocation length are shown in each figure. 
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Figure 6. Positions of the binding sites of H atom around the core for (a) easy core and (b) hard 
core configurations. Circles and triangles represent binding site with binding energy larger and smaller 
than 100 meV, respectively. Labels below the symbols denote the names of the binding sites. Electron 
density's minimum value along the screw dislocation direction is also shown by the color shades and 
contours. The high density peaks correspond to the atomic rows. 



